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ABSTRACT 

A new set of accurately measured frequencies of solar oscillations are used to infer 
the rotation rate inside the Sun, as a function of radial distance as well as latitude. 
We have adopted a regularized least squares technique with iterative refinement for 
both 1.5D inversion using the splitting coefficients and 2D inversion using individual m 
splittings. The inferred rotation rate agrees well with earlier estimates showing a shear 
layer just below the surface and another one around the base of the convection zone. 
The tachocline or the transition layer where the rotation rate changes from differential 
rotation in the convection zone to almost latitudinally independent rotation rate in 
the radiative interior is studied in detail. No compelling evidence for any latitudinal 
variation in position and width of tachocline is found though it appears that the 
tachocline probably shifts to slightly larger radial distance at higher latitudes and 
possibly becomes thicker also. However, these variations are within the estimated 
errors and more accurate data would be needed to make a definitive statement about 
latitudinal variations. 
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1 INTRODUCTION 

The measured splittings of solar oscillation frequencies of- 
fer us a valuable tool for studying the rotation rate in- 
side the Sun. It is possible to obtain both radial and lat- 
itudinal variation in the rotation rate (Schou, Christensen- 
Dalsgaard & Thompson 1994). Various techniques have been 
employed for inverting the splitting coefficients or even in- 
dividual splittings in a multiplet (Brown et al. 1989; Gough 
& Thompson 1991; Pijpers & Thompson 1992; Sekii 1993; 
Wilson & Burtonclay 1995; Corbard et al. 1997). 

The results obtained so far suggest that the observed 
surface differential rotation of the Sun persists through the 
convection zone (CZ). The rotation rate is nearly constant 
along different latitudes in most of the CZ, while in the 
radiative interior it is almost like rigid body rotation with 
a value intermediate between that of the solar equator and 
pole at the surface (cf., Thompson et al. 1996; Kosovichev et 
al. 1997). The transition occurs over a fairly thin layer, which 
is referred to as the "tachocline" (Spiegel & Zahn 1992). The 
thickness of the transition layer seems to be smaller than the 
best resolution that is currently achieved by inversion meth- 
ods. This layer contains a substantial radial gradient of ro- 
tation velocity, of opposite signs in low and high latitudes. 
It is widely believed that the tachocline with its angular 
velocity gradients could be the seat of the dynamo responsi- 



ble for the solar magnetic cycle (Weiss 1994; Gilman & Fox 
1997). The introduction of a toroidal magnetic field in this 
layer with latitudinal differential rotation is naturally ex- 
pected to lead to interesting consequences for the operation 
of the dynamo and its resultant manifestation in the solar 
surface activity. The strong gradient in the rotation rate is 
also expected to produce turbulence which is likely to mix 
material just below the convection zone — a phenomenon 
needed to match the structure of solar models with the he- 
lioseismically determined structure of the Sun (Richard et 
al. 1996; Basu 1997). The accurate measurement of solar in- 
ternal rotation rate also provides strong constraints on the 
theory of angular momentum transport in the stellar inte- 
rior (Rudiger & Kitchatinov 1996), which should contribute 
to our understanding of the solar spin down over its lifetime. 

In this work we investigate the internal rotation rate 
of the Sun, with particular emphasis on the tachocline. 
We have adopted the 1.5D inversion technique for invert- 
ing the rotation rate from the measured splitting coefficients 
(Ritzwoller & Lavely 1991; Schou, Christensen-Dalsgaard & 
Thompson 1994; Antia & Chitre 1996). We have also used a 
two-dimensional inversion of the individual frequency split- 
tings themselves. 

The location and structure of the tachocline is thus cru- 
cial in many models of the solar dynamo and it has been the 
subject of several detailed studies. Using a simple forward 
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modelling Kosovichev (1996) found that the tachocline is 
centered at a radial distance of (0.692 ± 0.005)i?o and with 
a width of (0.09 ± O.O4)i? , while Charbonneau et al. (1997) 
found it to be centered at (0.704 ± O.OO3)i?0 and with a 
width of (0.050 ± 0.012)_Rq. Similar results were also found 
by Basu (1997), who found that the tachocline is centered at 
(0.7050±0.0027)i? Q with a half-width of (0.0098±0.0026)i? Q 
which is (0.0480 ± O.O127)i? when scaled to the width as 
define by Kosovichev (1996) and Charbonneau et al. (1997) 
(see Section 3.1 for the definition of the width). Wilson, 
Burtonclay & Li (1996) find the tachocline to be some- 
what deeper at r = 0.68 ± 0.01-Rq and also slightly thicker 
(O.I2.R0). While all these results are roughly in agreement 
with one another, the differences in thickness are quite sig- 
nificant from the point of view of dynamo models as well as 
the hydrodynamical stability. Similarly, the exact location 
of tachocline with respect to the base of the convection zone 
is also crucial. Further, all these studies effectively assumed 
that the position and thickness of tachocline are indepen- 
dent of latitude. In this work, we attempt to find latitudinal 
variations in properties of the tachocline and also use im- 
proved data from GONG network to obtain better estimate 
for the tachocline properties. We use forward modelling tech- 
niques to detect possible latitudinal variation in properties 
of the tachocline. This includes the calibration method used 
by Basu (1997) and another technique based on simulated 
annealing. 

For this work we have used a number of different data 
sets obtained by the Global Oscillations Network Group 
(GONG) project (Hill et al. 1996). These data are very pre- 
cise and scan a large range of frequency and degree of modes. 
Apart from this we also use the data from BBSO (Woodard 
& Libbrecht 1993) combined with splitting coefficients for 
low degree modes as measured by BiSON (Elsworth et 
al. 1995). 

The rest of the paper is organized as follows. In Sec- 
tion 2 we outline the methods used to invert the data to ob- 
tain the solar rotation rate and describe the inversion results. 
Techniques for determining whether there is any latitudinal 
variation in the tachocline are summarized in Section 3. In 
Section 4 we discuss results of inversion when the contri- 
bution from the tachocline is removed from the data before 
inversion and our conclusions are stated in Section 5. 



Sun by 



2 INVERSIONS TO DETERMINE THE 
ROTATION RATE 

The different modes of solar oscillations can be described by 
three integers: the radial order n, the angular degree i and 
the azimuthal order m. The integers i and m are the de- 
gree and order respectively of the spherical harmonic func- 
tion used to describe the angular behaviour of the mode. 
In a spherically symmetric, non-rotating star, the frequency 
w n ,e,m of an eigenmode is independent of m and the mode is 
(21 + l)-fold degenerate. The spherical symmetry of the Sun 
is broken by rotation, lifting the degeneracy of the modes. 
The differences in frequency of modes with the same n and 
£, but different m, can be related to the rotation rate in the 



&n,e,m — ^>n,t,—m 

2m 

/ drdcos9 K nAm (r,9)Q(r,9), 
J -1 

where the kernels K n ^ im (r,8) are defined by Pijpers (1997). 

Most helioseismic data sets do not contain frequencies 
of individual modes or the individual splittings D n ^ >m as 
defined by equation (1), but rather frequencies of modes for 
a given (n,£) are expressed as sum of polynomials in m, 
namely, 



0J n ,t,m — ^n,l + 



E 



(2) 



where V^p (m) are suitable polynomials of degree s, and gen- 
erally, s max < 21. For a proper choice of the polynomials, 
the individual inversion problems for each splitting coeffi- 
cient Cs™'^ becomes decoupled from the rest (Ritzwoller & 
Lavely 1991). 

The data from the GONG instrument are available as 
frequencies of all individual modes, as well as splitting co- 
efficients for the polynomials as defined by Ritzwoller & 
Lavely (1991). We adopt both these forms of data, and use 
the so-called 1.5D inversion method (described below) to 
invert the data in the form of splitting coefficients. This 
method has the advantage of being efficient in terms of com- 
puting resources. However, in order to exploit the full poten- 
tial of the data we need to invert the individual frequency 
splittings directly using a two dimensional inversion method. 

2.1 The 1.5D inversion 

The rotational splitting coefficients are sensitive only to the 
component of rotation velocity that is symmetric about the 
equator and we therefore assume the rotation velocity to be 
symmetric. In order to determine the latitudinal dependence 
of the rotation rate, we follow Ritzwoller and Lavely (1991) 
and express the rotation velocity as 



v rot (r,9) = Sl(r,0)rsmO-- 



00 ~ 
= -5> 2s+ i(r)-y 2 ° s+ i(0), 

s=0 



(3) 



where 9 is the colatitude, Y£ (9) are the spherical harmonics 
and w s (r) are expansion coefficients which are related to the 
splitting coefficients cf" (cf., equation (2)) by 



S CM) = f 
Jo 



w s (r)K { ™' l) (r)r 2 dr, 



(4) 



where the kernels Ki"' £ \r) are given by (Ritzwoller & 
Lavely 1991) 



K 



(nJ) 



P0 



ie + i(i + ml - 2£r& - H s + ml) 



Jo Po(er+£(i+m> 2 dr 



(5) 



Here, po(r) is the density in the equilibrium solar model, 
while £ r and £/j are respectively, the radial and horizontal 
components of displacement eigenfunctions. Using the split- 



ting coefficients Cs ' from the GONG data, equation (4) 
can be inverted to obtain w s (r). The advantage of this choice 
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for expansion is that the resulting inverse problems for de- 
termining the individual components wi(r), ws(r),... get 
decoupled and each component can be estimated indepen- 
dently. The components w s (r) are calculated by solving sep- 
arate one- dimensional inversion problems with the iterative 
refinement of the regularized least squares solution (Antia, 
Chitre & Thompson 1996). The rotation rate at any given 
radial distance and colatitude can then be computed using 
equation (3). 

We use cubic B-spline basis functions to represent the 
rotation rate and the regularized least squares inversion is 
performed using the singular value decomposition. The B- 
splines are defined over a set of 50 knots which are uniformly 
spaced in acoustic depth. We have used only the first 6 terms 
of the expansion (3), as the higher splitting coefficients in the 
GONG data appear to be dominated by random noise. The 
observed rotational splitting coefficients from GONG data 
for t < 150 and 1 < v < 3.5 mHz are used for inversion. 
Further, in actual practice we directly find the individual 
components of rotation rate as defined by 



n s (r) 



2s + 1 w s (r) 
4-7T r 



(6) 



instead of w s (r). 

Since the inversion problem defined by equation (4) is in 
general ill-conditioned, some regularization or smoothing is 
required to obtain any meaningful solution in the presence of 
errors in observed data sets. In order to study the sensitivity 
of inversion results on smoothing prescription, we have tried 
three different prescriptions for the regularized least squares 
inversion: (i) First derivative smoothing, where we minimize 



K^™'^ (r)w s (r)r 2 dr 



+ A 



Jo r\dr) 



dr, 



(7) 



Here ci"' f ' are the observed splitting coefficients and ai n,e ^ 
the corresponding error and A is the regularization parame- 
ter. For this case, the rotation rate in the solar core — where 
the amount of information is rather meager as the splitting 
data for low degree modes have large errors — tends to a 
constant value when sufficient smoothing is applied. 

(ii) Second derivative smoothing, where we use the sec- 
ond derivative instead of the first in the second term of equa- 
tion (7). In this case, the rotation rate in the solar core tends 
to a monotonic linear profile, which is perhaps unrealistic 
as it may keep rising or falling depending on the gradient. 
In order to overcome this problem we apply the following 
boundary conditions at the center: 



0, 



d£h 
dr 

f2i(0) = 



(i > 1). 



(8) 



The second boundary condition may have some justification, 
as there is no information available to determine the higher 
order coefficients Qi(r) (i > 1) in the solar core from the 
splitting data. This is because only the first splitting coeffi- 
cient ci is known from observations for the low degree modes 
which sample the solar core . This condition is also required 
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Figure 1. The first three components of the rotation rate ob- 
tained with different prescriptions of smoothing. In all three pan- 
els the continuous line is the result obtained with first derivative 
smoothing. The short-dashed, long-dashed and dot dashed line 
are for second-derivative smoothing with the boundary conditions 
(equation (8)) applied at r = 0, 0.1 and O.3i?0 respectively. The 
dotted lines show la errors on solution obtained with the second 
derivative smoothing and boundary conditions applied at 0.3-Rq. 

to ensure regularity of rotation velocity at the origin (Cor- 
bard et al. 1997). 

(iii) This is same as (ii) except that the boundary con- 
ditions given by equation (8) are applied at a radial distance 
r = 0.3.R© (or r = O.IRq). These boundary conditions es- 
sentially try to constrain the rotation rate to become con- 
stant in the region inside where the boundary conditions are 
applied. This may be justified as we do not appear to have 
enough information to determine the gradient in rotation 
rate in the core (cf., Chaplin et al. 1996). 

In all these cases, the method of iterative refinement ef- 
fectively chooses the regularization parameter A as explained 
by Antia, Chitre & Thompson (1996). It is found that the 
regularization parameter increases with the order s of the 
splitting coefficient, which probably reflects the fact that 
higher order coefficients are dominated by noise in most re- 
gions. 



2.1.1 Results of 1.5D inversion 

We use the 1.5D inversion technique outlined above to in- 
fer the rotation rate in the solar interior using the GONG 
data for months 4-14, which consists of splitting coefficients 
for modes with I < 150. The \ 2 P er degree of freedom is 
found to be close to unity (between 1.1-1.2) in all cases. In 
order to study the influence of smoothing on the inversion 
results we perform inversion using different smoothing pre- 
scriptions given in Section 2.1 and the results are displayed 
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in Fig. 1. This shows the rotation rate corresponding to the 
first three splitting coefficients ci, C3 and c 5 . From the fig- 
ure it is clear that the results in the convection zone are 
not very sensitive to the choice of smoothing or to the point 
at which the boundary conditions given in equation (8) are 
applied. Noticeable differences are seen only for those parts 
of the Sun where the data have large errors. The maximum 
difference of about fO nHz occurs in the core in the latitu- 
dinally independent component fii. Similar differences are 
seen for the component ^3 also. These differences are com- 
parable to the error estimate in the inversion results arising 
from errors in observed splitting coefficients. It may be noted 
that the error estimates shown in the figure for the second 
derivative smoothing with the boundary conditions applied 
at r = 0.3i?© show a decrease in the core. This is artificial, 
and is entirely the result of boundary conditions. In reality 
the errors should increase rapidly as r decreases in the core. 
Some of the differences in the core between the results using 
first and second derivative smoothing are due to absence of 
boundary conditions in the first derivative case. However, a 
part of the difference may also be due to possible system- 
atic errors in splitting coefficients of the low degree modes. 
From this figure it is clear that for r > 0.5Rq the choice of 
smoothing or the point where the boundary conditions are 
applied does not make significant difference and in most of 
the work we have confined ourselves to this region. All the 
following results using 1.5D inversion have been obtained 
using second derivative smoothing with the boundary con- 
ditions applied at 0.3Rq. 

In order to study the sensitivity of the inversion tech- 
nique to possible systematic errors in the input data sets 
we have repeated the inversion for various sets of GONG 
data and the results are shown in Fig. 2. This figure also 
includes the results obtained using the averaged BBSO data 
for years 1986, 1988, 1989 and 1990 (Woodard & Libbrecht 
1993), combined with the splittings for low degree modes 
from BiSON (Elsworth et al. 1995). It is clear that there is 
some systematic difference between different data sets, but 
the difference is comparable to the error estimates arising 
from inversions. This difference is also comparable to that 
arising from different smoothing prescriptions as displayed 
in Fig. 1. The spread between various curves in Fig. 2 should 
give an estimate of expected errors in inverted profiles, in- 
cluding those arising from systematic errors in the input 
data. It may be noted that the data from GONG month 10, 
GONG months 4-7 and BBSO+BiSON have larger errors as 
compared to that in the GONG months 4-14 data and this 
is reflected in tachocline region where the GONG months 
4-14 data appear to have higher resolution and hence the 
tachocline is sharper. 

A contour diagram showing the rotation rate inside the 
Sun as inferred using the GONG months 4-14 data is shown 
in Fig. 3. It can be seen that the rotation rate is approxi- 
mately constant along radial lines in the convection zone. 
These results are similar to earlier inversions for rotation 
rate (Thompson et al. 1996; Kosovichev et al. 1997). The 
tachocline is clearly visible in these contour diagram. Apart 
from the tachocline there, is another shear layer near the 
solar surface where the rotation rate increases with depth. 
This shear layer appears to extend to all latitudes and the 
rotation rate increases by about 17 nHz in this layer at the 
equator, but the change is lower at higher latitudes. The 
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Figure 2. Solar rotation rate at the equator, 30° and 60° lati- 
tudes obtained using different data sets. In all three panels, the 
continuous lines arc the results obtained using the GONG months 
4—14 data with the dotted lines showing the 1<t error limits. 
The short-dashed, long-dashed and dot-short dashed lines are for 
GONG months 4-10, months 4-7 and month 10 data respectively, 
while the dot-long dashed line is for the BBSO+BiSON data com- 
bination. Note that while for all the GONG data sets we have used 
splitting coefficient from ci to en, for the BBSO+BiSON set we 
have used only the data up to C5. 

maximum value of rotation rate occurs around r = O.95i?0. 
The maximum rotation rate at the equator is 467.5 ± 0.2 
nHz at r = 0.9457?©. The inverted rotation rate at the solar 
surface is close to that inferred from Doppler measurements 
(Snodgrass 1992). The rotation rate in the radiative interior 
is more or less constant and some of the features seen in the 
contour diagram do not appear to be significant. Although 
there is considerable uncertainty in the estimate of the ro- 
tation rate in the core, but it appears to be less than the 
surface equatorial rotation rate. 

2.2 2D inversion 

Although the 1.5D inversion technique described in Sec- 
tion 2.1 is very efficient in terms of computing resources, 
it is not clear if the expansion of rotation rate given by 
equation (3), imposes any limitation on the solution. A pos- 
sible drawback of 1.5D inversion is the loss of information, 
since the number of splitting coefficients is generally much 
smaller than the number of individual splittings, but it is not 
clear if the individual splittings contain any more informa- 
tion than the first few splitting coefficients. A more serious 
problem occurs in the inversion of higher order coefficients, 
which have useful information only in the outer part of the 
convection zone and as a result the solution in most of the 
interior is essentially determined by the applied smoothing 
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Figure 3. A contour diagram of the solar rotation rate as ob- 
tained by the 1.5D inversion technique using GONG months 4-14 
data. Due to the symmetry of the inversion results, the rotation 
rate has been shown for just one quadrant only. The dotted con- 
tours have been drawn at intervals of 5 nHz., and the continuous 
ones at intervals of 20 nHz. The thick continuous line is the con- 
tour at a level of 440 nHz. The x-axis represents the solar equator 
while the j/-axis represents the rotation axis. 



and boundary conditions. A 2D representation of rotation 
rate will hopefully be able to overcome this problem. An- 
other drawback of the 1.5D inversion is that all components 
£l s (r) make their maximum contribution at the pole and fur- 
ther this maximum value is much larger than that in other 
regions. This is particularly true for the higher order com- 
ponents. As a result, the errors in the 1.5D inversions tend 
to get highly magnified near the pole and may even give rise 
to spurious features if very high order terms are included. In 
order to overcome these problems we attempt a 2D inversion 
technique which does not use the expansion given by equa- 
tion (3) , but instead directly represents the two dimensional 
function fl(r, 9) in terms of suitable basis functions. 

In order to solve the inversion problem defined by equa- 
tion (1) we represent the rotation rate in terms of B-spline 
basis functions in r and 9, 



n(r,0) = ]T2 6 y* i ( r )^( COB *)' 



(9) 



= 1 3 = 1 



where bij are the coefficients of expansion and <j>i(r) are the 
B-spline basis functions over r and tpj(cos6) are those over 
cos#, n r and n$ are the number of basis functions in r and 
cos 9 respectively. We use a set of knots which are uniformly 
spaced in acoustic depth and cos9 respectively to define 
4>i(r) and i/)j(cosS). 

This inversion problem is also solved using the regu- 
larized least squares technique with iterative refinement. In 
this case, we have used only second derivative smoothing, 




Figure 4. A contour diagram of the solar rotation rate as ob- 
tained by the 2D inversion of individual splittings using GONG 
months 4-14 data. The format is the same as that for Fig. 3. 

which involves minimizing 



n,£,m 



+ A, 
+ A £ 



/■«© r L 

D ni t }m - / dr d cos 6K n ,t tm {r,0)£l{r, 6) 

Jo J -l 



d cos 9r 



d 2 Q. 
dr 2 



d 2 n 

d cos 9 2 



where, A r and Xe are the two regularization parameters con- 
trolling the smoothing. No boundary conditions are applied 
in this case to constrain the rotation rate in the core. We 
have used 50 knots in r and 30 knots in cos# to represent 
the rotation rate. 

It is also possible to perform 2D inversion for split- 
ting coefficients, (Schou, Christensen-Dalsgaard & Thomp- 
son 1994; Pijpers 1997) where the rotation rate is expressed 
in terms of 2D basis functions (equation (9)) and appro- 
priate combinations of individual splittings are constructed 
to relate the corresponding splitting coefficients to the ro- 
tation rate. In order to see whether the differences between 
the 1.5D and 2D results are due to the expansion of rotation 
rate or the data, we have done a two-dimensional inversion 
for the splitting coefficients also. Thus we have two sets of 
results using 2D inversions, one for 2D inversion of individ- 
ual splittings D n> ( im , and the other for 2D inversion of the 
splitting coefficients ci"' e \ 

2.2.1 Results of 2D inversion 

A contour diagram of the rotation rate inferred by 2D inver- 
sion of GONG months 4-14 data, comprising of about 85000 
splittings of individual modes is shown in Fig. 4. The x P er 
degree of freedom in this case is around 1.3. Note that inside 
the CZ, the results are essentially similar to that obtained 
by the 1.5D method, despite having a much larger number of 
splittings in 2D inversion. Thus it appears that the data are 
well represented by the 6 splitting coefficients. However, in 
the radiative interior the solutions obtained using 1.5D and 
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Figure 5. Comparison of the rotation inversion results at fixed 
latitude obtained by the 1.5D and 2D inversion methods. The 
continuous line shows the results obtained by the 1.5D inversion 
method with the dotted lines showing the la error limits. The 
dashed line are those obtained by the 2D inversion of the individ- 
ual splittings and the dot-dashed line are those obtained by the 
2D inversion of the splitting coefficients. 

2D inversions are significantly different. A large part of the 
difference is due to the boundary conditions imposed in the 
1.5D inversion which along with the smoothing, tend to pro- 
duce solid body rotation in the interior. In the absence of any 
boundary condition, the 2D inversion technique attempts to 
fit the splittings for low degree modes which probably have 
some systematic errors, and produces a sharply decreasing 
rotation rate in the radiative interior. This decrease may not 
be real as it could result from second derivative smoothing 
coupled with errors in data. 

It appears that the behaviour of the solutions in the po- 
lar regions is also somewhat different. The 2D-solution shows 
a rapid decrease of the rotation rate towards the pole similar 
to that seen in the data from the MDI instrument on board 
the SOHO spacecraft. The 1.5D result also shows a decrease 
in the surface rotation rate at the pole, but the reduction is 
not as much as in 2D inversion of individual splittings. The 
errors in inversion increase rapidly with latitude near the 
pole and as a result it is difficult to discern any features at 
high latitudes reliably. Thus the differences at high latitudes 
possibly reflect our inability to obtain reliable inversion re- 
sults in that region. However, from the form of expansion of 
rotation rate (equation (3)) in 1.5D inversion, it is clear that 
contribution of each component (fij) has a strong maximum 
at the pole and any error in these components will be highly 
magnified there. This is particularly true of the higher order 
terms in the expansion. It thus appears that 2D inversion of 
individual splittings, which does not assume any particular 
expansion of rotation rate, may be able to give better re- 



c 
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Figure 6. Comparison of the rotation inversion results at fixed 
radius. The line styles are the same as in Fig. 5. In the top panel 
the heavy short dashed- long dashed line shows the observed sur- 
face rotation rate as estimated by Dopplcr measurements (Snod- 
grass 1992). 

suits near the pole, though the errors will still be large and 
the smoothing will play a dominant role in determining the 
solution near the poles. 

In order to see whether the differences between the 1.5D 
and 2D results are due to the representation of rotation rate 
or the data, we have also done a two dimensional inversion of 
the splitting coefficients. For this we use the same splitting 
coefficients as were used in the 1.5D inversion but expand the 
rotation rate in terms of 2 dimensional basis functions using 
equation (9). These results are also shown in Figs. 5 and 6. 
Note that in the outer layers of the Sun the results of the 
1.5D and two-dimensional inversion of splitting coefficients 
are very close to each other, but slightly different from the 
results of the 2D inversion of the individual splittings. In the 
deeper layers however, the results of 1.5D inversion are quite 
different from both the 2D inversions. As already discussed 
this is due to the boundary conditions (equation (8)) that 
are applied in the 1.5D inversion method. 



3 THE TACHOCLINE 

It can be seen from the results of the inversions that there 
is a sharp transition close to the base of the convection zone 
where the rotation rate changes from differential rotation 
in the convection zone to a rotation rate that is almost in- 
dependent of latitude. Since the inversions are not able to 
resolve the tachocline, other techniques have been employed 
to study this shear layer. Even though inversions indicate 
that the tachocline is slightly shallower and thicker at high 
latitudes, it is not clear if this represents a real variation or 
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Figure 7. The splitting coefficients for GONG months 4—14 com- 
bined to obtain the data for different latitudes, which arc marked 
in the figure. The points represent the combinations binned in 
groups of 15. 

is just an artifact of inversion technique caused by the fact 
that the extent of jump increases with latitude and the res- 
olution of inversion techniques deteriorates with increasing 
latitude. For an independent confirmation of this variation 
we construct combinations of rotational splitting coefficients 
ci n ' e \ which give the rotation velocity at some predefined co- 
latitude #o- Thus multiplying equation (4) by dY® / d8\g = g 
and summing over s we get 



dYl. 



d6 



/■«© 

/ Vr, 

Jo 



■, ot (r,e )K {n ' l) {r)r z dr.(ll) 



Here, in principle, the kernel K^ t,n \r) also depends on s, 
but for simplicity we neglect the s dependence, since the 
s-dependent term is in general very small, as can be seen 
from equation (5). With this approximation equation (11) 
reduces to an one dimensional inversion problem at fixed 
latitude. Since inversion cannot resolve the tachocline, we 
use a forward modeling technique to estimate the parameters 
defining the tachocline. 

We have first made the appropriate combinations of the 
splittings data for different latitudes using the coefficients 
ci-cn. Fig. 7 shows the data, binned in groups of 15 modes, 
plotted as a function of the lower turning point. We show 
the data at the equator, 45° and 60° latitudes. Note that 
the equator shows some hint of a jump, while 45° and 60° 
latitudes show a clear jump. We concentrate on the latitudes 
which show evidence of the transition and try to determine 
the magnitude of the jump, position and thickness of the 
transition at latitudes of 0°, 15° , 45°, 60° and 75°. The data 
for the 30° latitude does not give any clear indication of the 
transition, which is consistent with inversion results. The 



latitudinal variation of the magnitude of the jump is obvi- 
ous from the figure, but the change in position and thickness 
are not clear. We use both the calibration method used by 
Basu (1997) and another method based on simulated an- 
nealing to study the tachocline at each latitude separately. 

Apart from this we have also tried a 2D fit with some 
assumed form for tachocline with latitudinal variation to 
simultaneously fit all the splitting coefficients for obtaining 
the latitudinal variation in the properties of the tachocline. 

3.1 Calibrating the tachocline 

The method followed here is the same as that used by Basu 
(1997). Data on frequency splittings for modes trapped in 
the convection zone and those with turning points around 
the CZ base are quite precise. The exact position and thick- 
ness of the tachocline can be determined by calibrating 
the difference in splittings between the Sun and models 
with known position and thickness of the tachocline. While 
Basu (1997) assumed that the position and thickness of the 
tachocline can be determined by the splitting coefficient C3 
alone, which automatically ensures that the position and 
thickness of tachocline are independent of latitude. In this 
work, we explicitly study the latitudinal variation by ap- 
plying the same procedure to data for each of the chosen 
latitudes. 

In order to estimate the position, thickness and the 
jump in the tachocline we construct a series of models with 
known properties. We parameterize the calibration models 
with the rotation profile: 

^cal = — ^ ry— r, (12) 

1 + exp[(r d - r)/w\ 

where 8Q is the jump in the tachocline, w is the half-width 
of the transition layer, and r d the mid-point of the transition 
region. Thus the rotation rate increases from a factor 1/(1 + 
e) of its maximum value to the factor 1 — 1/(1 + e) of its 
maximum value in the range r — r d — w to r — r d + w. The 
models described by equation (12) have almost zero rotation 
rate in the core, which is, of course, not true for the Sun. 
To take this into account we subtract the contribution of 
a uniform rotation rate fi c , the estimated value of rotation 
rate in the interior as obtained from the inversion results, 
from the observed splittings. 

Note that our parameterisation of the tachocline is dif- 
ferent from that of Kosovichev (1996) and Charbonneau et 
al. (1997). The definition of the position and the jump re- 
mains the same, the thickness of the tachocline as defined 
by them is roughly 4.9 times the half-width we have defined, 
i.e. a half-width of O.OITi© in our model corresponds to a 
thickness of O.O49i?0 in their models. Thus the tachocline 
thickness of (0.050 ± O.O12)i?0 as estimated by Charbon- 
neau et al. (1997) using their model will be equivalent to a 
half-width of (0.0102±0.0025))i? Q by our definition, which is 
consistent with the value of (O.OO98±O.OO26)J?0 determined 
by Basu (1997). Of course, the form of variation in rotation 
rate inside the tachocline can only be verified by inversions, 
which do not at present have the required resolution. How- 
ever, as long as the thickness of tachocline is small enough, 
the exact form of variation within this layer may not be im- 
portant, as is seen by the similarity of the results obtained 
by Charbonneau et al. (1997) and Basu (1997). 
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If ce be the combination of observed splitting coeffi- 
cients for a given latitude (after removing the contribution 
from O c ) and ae that of the calibration model, then to de- 
termine the jump, we make the following fit by treating the 
splittings as a function of the lower turning point, r t , of the 
mode: 

cg(r t ) = aa g (r t ) + 4>(r t ), (13) 

where a is the factor determining the jump and 4>(rt) is a 
low degree polynomial which takes into account any trend 
in the real rotation rate not taken into account by the pa- 
rameterisation in equation (12). The constant term in this 
polynomial will also take care of differences arising due to 
use of incorrect £l c while subtracting out the contribution 
from core rotation rate. A polynomial of degree 2 is found 
to be sufficient. The fit is made between r t of 0.6 and 0.9 
Rq . Modes with lower turning point are avoided because of 
large observational errors in these modes. Modes with higher 
rt are not used so that shear layer known to exist just below 
the solar surface does not affect the results. We perform a 
least squares fit with the weights for each mode being the 
inverse of the errors in the corresponding splitting. 

We follow exactly the same procedure as that of Basu 
(1997) for determining the position and thickness of the 
tachocline. To recapitulate briefly, the difference in the split- 
ting coefficients between a model and similar models which 
have discontinuities at different positions have a well defined 
peak and the height of the peak is proportional to the dif- 
ference in the positions of the discontinuity. Thus the peak 
height can be calibrated to find the position of the disconti- 
nuity. If the models are not similar, e.g., have different trends 
in the convection zone or in the interior, or have a different 
width of transition, the peak lies on a smooth part which 
can be represented as a low degree polynomial. A similar 
calibration can be used to estimate the thickness over which 
the transition of the rotation rate occurs. In this case the 
width of the peak between two models is proportional to the 
thickness of the wider transition, but a simple scaling of the 
radius around the peak position can reduce the curves to 
similar widths. 

We have constructed models with r<j of 0.68, 0.69, 0.70, 
0.71 and 0.72 R Q and half-width w of 0, 0.005, 0.01, 0.015, 
0.02 and 0.025 Rq. We have used SSI of 20 nHz in the cal- 
ibration models and determined the actual jump from the 
splittings as outlined earlier. The splittings in the calibration 
models are then scaled to the required jump. 

For the purpose of calibration, we consider the differ- 
ence in the splitting coefficients between neighbouring cali- 
bration models and fit a spline through the points 

<t>{r)=5a{r) = ^a i il> i {r), (14) 

i 

where the $(r)'s are the calibration curves and ip(r) are the 
cubic B-spline basis-functions. Thus we have 4 calibration 
curves from the 5 calibration models. 

The difference in splittings between each calibration 
model and the observed splittings can be fitted with the 
form 

Sc e = cs(r) - Ofl(r) = a$(r) + f(r). (15) 

Here 3>(r) is the calibration curve defined in equation (14) 
and f(r) is a low degree polynomial used to take into account 



systematic effects arising from differences in other parame- 
ters, like width etc. between the observations. As in Basu 
(1997) we find that a polynomial of degree two or three is 
optimum. The constant a and the coefficients of the poly- 
nomial f(r) are obtained by a least-squares fit to the data. 
In practice, we determine a for all five calibration models 
and interpolate to find the points where a = 0. The four 
calibration curves give four results which are then averaged. 

Ideally, SQ, ra, and w should be determined simulta- 
neously, however, for simplicity in this work we determine 
these parameters by independent fits. It has been shown by 
Basu (1997) that statistical errors arising from uncertainties 
in observed splittings dominate over the systematic errors, 
and therefore, such a procedure may not introduce signifi- 
cant additional errors. To try and keep the parameters of the 
calibration models close to that of the actual tachocline, the 
fit in practice is done in two steps. We first determined the 
positions and widths of the tachocline at the different lati- 
tudes using the models with parameters found for the coeffi- 
cient C3 by Basu (1997). The process was repeated with cal- 
ibration models constructed with parameters closer to those 
determined in the first round of fits. 

Note that the models defined by equation (12) have 
a fiat rotation rate in the CZ, which is obviously not the 
case at all latitudes. This is taken care of in our fitting 
process described above through the smooth part, <f>(r) in 
equation (13) and f(r) in equation (15). However, in order 
to check how much difference the flat rotation rate in the CZ 
makes, we have also constructed calibration models where 
the rotation rate in the CZ follows the trend revealed by the 
inversions described in the previous section. In addition to 
the form shown in equation (12), these models have latitude 
dependent extra terms, with rotation rate defined as: 

of ftcai + B(r - 0.7) if r < 0.95 , . 

\ fi ca i - C(r - 0.95) + 0.25B if r > 0.95 1 ' 

Here, the coefficients B and C are obtained from the inver- 
sion results and the term f2 ca j is defined by equation (12). 
The models described by equation (16) use an approximate 
value of the jump; we still use the fit in equation (13) to find 
the exact magnitude of the jump. The procedure to find the 
position and thickness of the tachocline remains the same. 
It must be noted that for these models, the polynomial <j) of 
equation (13) and f(r) in equation (15) are very small. 

For all cases, the error estimates in the tachocline pa- 
rameters are obtained using Monte-Carlo simulations. We 
have used two sets of GONG data, those from GONG 
months 4-10 and GONG months 4-14 for the present work. 



3.2 The method of simulated annealing 

The calibration method described above suffers form the dis- 
advantage that each parameter defining the tachocline is de- 
termined separately when the rest of the parameters are held 
fixed. We therefore tried another forward modeling method, 
where the jump, position and width can be found simultane- 
ously. Since this will require a nonlinear least squares fit, we 
have resorted to the method of simulated annealing, which 
has better chance for finding the global minimum. For this 
purpose the rotation rate at any given latitude is parame- 
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terized by 



^ann (f) 



n c + B(r - 0.7) 



_l m 

l+exp[(r d -r)/tu] 

tic + 0.25B - C(r - 0.95) 



+ 



l + cxp[(r d -r)/u>] 



if r < 0.95 



if r > 0.95 



(17) 



where fi c , B, C are the three parameters defining the 
smooth part of rotation rate while <5f2, rd and w define the 
tachocline. Here B is the average gradient in the lower part 
of convection zone, while C is the gradient in the near sur- 
face shear layer. These six parameters are determined by a 
non-linear least squares fit to the combinations of splitting 
coefficients representing the rotation rate at the required 
latitude. Once again we use only those modes which have 
turning points in the range 0.6-0.97?© in this fit. We use the 
method of simulated annealing (Vanderbilt & Louie 1984; 
Press et al. 1993) to minimize the \ 2 function. Since there 
are likely to be many local minima where the minimiza- 
tion tends to get trapped, even with simulated annealing, 
we make 20 attempts using different sequence of random 
numbers in the annealing procedure to find the minimum 
and accept the one that gives the lowest \ ■ 

Instead of fitting the rotation rate at each latitude sepa- 
rately and then finding the variation in the tachocline prop- 
erties, we can directly fit a 2D form of rotation rate with 
tachocline to obtain this variation. The form fitted is the 
same as in the ID case (equation (17)) with the following 
substitutions 



B=B 1 + B 3 P i (6) + B 5 P 5 (6), 
8fl =8Q! + 6n 3 Ps{0) + 5n 5 P 5 (8), 
Td =Tdi + r d :iPi(6), 

W =Wi + W 3 Pz{6), 

where 

P 3 (6) = 5 cos 2 (9- 1, 

P 5 {6) = 21 cos 4 (9- 14 cos 2 6 + 1, 



(18) 



(19) 



are polynomials used to define the latitude dependence. We 
have used these polynomials so as to ensure some degree of 
separation between contributions to various splitting coef- 
ficients. This introduces 5 parameters to define the smooth 
part and another 7 parameters to define the tachocline. It is 
not clear if all these parameters are required to explain the 
data and we have carried out experiments involving different 
combinations to determine which of these parameters are re- 
quired. Once again we use the simulated annealing technique 
to simultaneously fit the first 3 splitting coefficients C1-C5 for 
all modes with lower turning point in the range 0.6-0.9-Rq . 
The reason for using only the first 3 splitting coefficients is 
that with the assumed form for tachocline model given by 
equation (18), we do not expect to fit the higher coefficients 
properly. If these coefficients are to be fitted then additional 
parameters will need to be introduced in <5H and B in equa- 
tion (18). We have tried this also but we do not think that 
the fits with additional parameters and coefficients are any 
better than the ones considered here as the higher order pa- 
rameters turn out to be rather small. As a result, in this 
work we present the results obtained by fitting only the first 
three splitting coefficients. 




20 40 60 80 

Latitude (°) 

Figure 8. A summary of the tachocline results for the test model. 
Panels (a),(b) and (c) show the results of the jump, position and 
half-width respectively. In each panel the continuous line repre- 
sents the exact value. The circles show the results of the cali- 
bration method and the triangles are the results obtained by ID 
annealing. The symbols are displaced by ±0.5° about the true 
latitude for the sake of clarity. In panel (c) the actual width has 
been scaled down by a factor of 3.25 to account for the difference 
in form of variation within the tachocline. 



3.3 Results 

In order to test the procedure outlined above we con- 
structed a test model with prescribed rotation rate including 
a tachocline and applied the techniques for determining the 
properties. The rotation rate in the test model was chosen 
to be 

Q(r,0) = 430 + 20 sin(7rr/2) 

+ 10(1 - 4 cos 2 6 - cos 4 6)Fj(r), 

{—1 if r < td — w 

sin(0.5-7r(r — rd)/w) if ra — w < r < ra + w (20) 
+ 1 if r > r d + w 

r d = (0.69 + 0.02 cos 2 6)R @ , 

w = (0.005 + 0.02 cos 2 6)R@. 

The splitting coefficients computed for this model were per- 
turbed by adding random errors with the same distribution 
as that specified by quoted errors in the GONG months 4- 
14 data. The perturbed data were then used to infer the 
characteristics of the tachocline and the results using the 
calibration and ID annealing methods are shown in Fig. 8. 
Since this model has a different form of variation within the 
tachocline as compared to the models we are using for the 
fits, we do not expect the thickness of tachocline as deter- 
mined by our procedure to agree with the actual thickness. 
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Comparing the region in which the rotation rate varies from 
1/(1 + e) to 1 — 1/(1 + e) of the total jump, it appears that 
the effective width in the test model is about 3.25 times 
less than what is given by equation (20), when models with 
form given by equation (17) are used. Thus the estimated 
values should be compared with this scaled width as has 
been done in Fig. 8. Apart from the form of variation within 
the tachocline, the trend in the lower CZ in test model is 
also far from linear and hence may not be properly mod- 
elled by the tachocline model used in annealing fits. Note 
that all the results are roughly within the error bars of the 
exact values, even though the form of variation within the 
tachocline as well as the trend in the CZ in the test model 
arc different from those in the calibration models. If we use a 
test model with the same form for tachocline as used in the 
calibration models it is possible to obtain much better re- 
sults. This gives us confidence that we can indeed determine 
the parameters of the solar tachocline. 

The variation in the position and thickness of tachocline 
with latitude is not totally clear from the estimated param- 
eters, since these variations were chosen to be comparable 
to the error estimates as happens to be the case for observed 
splittings also. It appears that a variation of 0.027?q in po- 
sition of tachocline is barely at the limits of detection with 
the present data. The variation in thickness is not very ev- 
ident as the estimated thickness appears to be too small at 
all latitudes. 

Having tested our techniques on a test model we now 
apply the same procedure to the GONG data for the months 
4-10 and 4-14. Fig. 9 shows the process of determining the 
tachocline parameters for the solar equator. The results ob- 
tained using the calibration methods are listed in Table 1. 
A positive jump implies a rotation rate which is higher in 
the CZ than in the radiative interior. From Table 1 we 
note that the results for the two data sets are consistent 
with each other within the estimated errors. However, the 
GONG months 4-10 data have larger error as compared to 
the months 4-14 data and that is reflected in the larger 
errors in the tachocline parameters. We thus focus our at- 
tention on results for data from GONG months 4-14, which 
have also been used in the method of simulated annealing. 

The change in the tachocline jump as a function of lat- 
itude is very clear. This is not surprising since the change is 
large enough to be seen by normal inversions also. The result 
obtained seems to depend somewhat on the type of calibra- 
tion model used, although at each latitude they are consis- 
tent within errors. It also appears that the results are more 
sensitive to data errors when models with a flat rotation- 
rate in the CZ are used. This is perhaps not surprising as 
the jump at each latitude is not assumed to be known be- 
forehand, while for the models with trend, a first estimate 
of the jump is made from the inversion results and only a 
correction-factor is obtained by the fit in equation (13). The 
estimated errors in position and thickness at all latitudes 
are larger than the corresponding errors in mean values as 
estimated by Basu (1997) using only the splitting coeffi- 
cient C3. This is partly due to the fact that the error esti- 
mates in splittings for each latitude is larger than those in 
C3 alone. Further, a large part of the variation in tachocline 
is determined by C3, which shows a clear jump around the 
tachocline, and at the same time has very little variation in 
the CZ or in radiative interior thus making it much easier 
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Figure 9. (a) Calibration for determining the jump in the ro- 
tation rate in tachocline. The points are the splitting-coefficients 
combination for the equator with the contribution from Q c re- 
moved. The continuous line is the fit to the data, which can be 
decomposed into two parts — the dotted line which is the term 
<j>(rt) in equation (13) and the dashed line which is the calibration 
model scaled to the fitted jump (term aag(rt))- The data used are 
the GONG4— 14 splittings. The dashed and continuous lines more 
or less coincide in this case, (b) The spline representation of the 
difference between the data shown above and the the five calibra- 
tion models used to determine the tachocline position. The con- 
tinuous, dotted, small-dashed, long-dashed and dot-dashed lines 
arc difference with calibration models for = 0.68, 0.69, 0.70, 
0.71 and 0.72Rq respectively. All these calibration models have 
w = O.OO5-f?0 (c) The spline representation of the difference be- 
tween the data and the five calibration models used to deter- 
mine the tachocline width. The continuous, dotted, small-dashed, 
long-dashed and dot-dashed lines are difference with calibration 
models for w = 0.005, 0.01, 0.015, 0.02 and 0.025.R© respectively. 
These calibration models have = 0.69-Rq . 

to fit the tachocline parameters. 

There appears to be a slight variation in the radial po- 
sition of the tachocline with the tachocline moving outwards 
at higher latitudes. However, the change is not very signif- 
icant - being only about la between the equator and 60° 
latitude. But the results seem to indicate a nearly systematic 
shift. 

The question about variation in the thickness is less 
clear, however. In fact for all latitudes, it appears that the 
thickness is very small and comparable to the error esti- 
mates. Thus better data with reduced errors are required 
before the thickness can be determined reliably. With this 
method, at the moment we can only put an upper limit on 
the thickness of the tachocline at all latitudes. Basu (1997) 
had shown that the thickness measurements are somewhat 
sensitive to the calibration models used. Thus we should try 
to check these results against those obtained by the tech- 
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Table 1. Solar tachocline parameters as determined by calibration 



Calibration models with flat CZ Calibration models with trend in CZ 



Lat. 


Jump 


Position 


Half-width 


Jump 


Position 


Half-width 


(°) 


(nHz) 


(Re) 


(Re) 


(nHz) 


(Re) 


(Re) 








GONG months 4- 


-10 data 









21.66 ± 2.57 


0.6921 ±0.0142 


0.0134 ±0.0120 


20.92 ± 2.46 


0.6954 ±0.0113 


0.0029 ± 0.0120 


15 


18.52 ± 1.91 


0.6941 ±0.0111 


0.0124 ±0.0137 


18.11 ± 0.92 


0.6891 ±0.0115 


0.0190 ± 0.0131 


45 


-27.81 ± 2.30 


0.6863 ± 0.0277 


0.0215 ±0.0173 - 


-27.17 ±2.06 


0.6903 ± 0.0263 


0.0170 ± 0.0172 


60 


-60.47 ±3.52 


0.7025 ±0.0117 


0.0070 ± 0.0075 - 


-60.08 ± 1.33 


0.6990 ± 0.0099 


0.0050 ± 0.0077 


75 


-86.41 ±9.99 


0.6968 ± 0.0224 


0.0094 ±0.0129 - 


-89.04 ±2.56 


0.6974 ± 0.0238 


0.0123 ± 0.0135 








GONG months 4- 


-14 data 









19.39 ± 2.00 


0.6944 ± 0.0096 


0.0079 ±0.0130 


21.52 ± 0.82 


0.6851 ±0.0077 


0.0047 ± 0.0083 


15 


18.13 ± 1.55 


0.6996 ± 0.0091 


0.0083 ± 0.0106 


18.29 ± 0.89 


0.6922 ± 0.0097 


0.0043 ± 0.0087 


45 


-28.89 ±2.66 


0.7077 ±0.0137 


0.0047 ± 0.0061 - 


-29.87 ±2.22 


0.7048 ± 0.0148 


0.0059 ± 0.0067 


60 


-57.18 ±4.11 


0.7058 ± 0.0098 


0.0031 ± 0.0055 - 


-57.10 ± 1.51 


0.7082 ± 0.0072 


0.0051 ± 0.0062 


75 


-88.21 ±6.70 


0.7204 ±0.0183 


0.0154 ± 0.0235 - 


-87.15 ± 1.78 


0.7162 ±0.0178 


0.0141 ± 0.0136 



Table 2. Tachocline parameters from ID-annealing 



Lat. 

(°) 



Jump 
(nHz) 



Position 

(Re) 



Half-width 

(Re) 




15 
30 
45 
60 
75 



17.20 ±4.96 
14.84 ±1.52 
-7.46 ± 1.71 
-33.83 ± 3.90 
-59.86 ± 6.19 
-91.19 ± 3.36 



0.6843 ±0.0112 
0.7146 ±0.0050 
0.7193 ±0.0196 
0.7160 ±0.0072 
0.7045 ± 0.0061 
0.6878 ±0.0089 



0.0020 ± 0.0019 
0.0028 ± 0.0022 
0.0242 ± 0.0136 
0.0122 ± 0.0073 
0.0089 ± 0.0078 
0.0216 ± 0.0095 



1.0401 
1.1271 
0.9860 
0.9276 
1.0712 
1.1762 



nique of simulated annealing. 

In order to obtain an independent measure of variation 
in properties of tachocline with latitude, we adopt the tech- 
nique of simulated annealing to fit the tachocline parame- 
ters to the GONG months 4-14 data for different latitudes. 
Fig. 10 shows the ID annealing result for the equator. Note 
that we get a good fit and the residuals arc random and con- 
sistent with the error estimates. The ID annealing results 
for various latitudes are listed in Table 2. As in the case 
of the calibration technique, the jump shows a clear change 
with latitude, although there appears to be some systematic 
difference between the value of the jump obtained by the two 
techniques. All the values appear to be reduced in annealing 
results as compared to the corresponding values obtained by 
the calibration method, although for individual latitudes the 
results are generally within error limits from those obtained 
using calibration method. The reason for this discrepancy 
is not altogether obvious, but there may be some ambiguity 
in the definition of jump as a part of the variation across 
the tachocline may be accounted for by the smooth trend, 
defined by the term involving B in equation (17). Note from 
the last column which gives the \ 2 P er degree of freedom, it 
is clear that the fit is reasonably good as all the values are 
close to unity. 

The results at 30° latitude are not particularly reliable 
as the jump is too small to define the tachocline properly and 
hence the errors are very large. At high latitudes although 
the splittings have larger error, the increase in magnitude 
of jump makes the tachocline better defined and in some 
cases the error estimate also reduces since the fits are more 
stable. At low latitudes, however, there is some problem in 
finding a proper fit, as the estimated values have a larger 
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Figure 10. The ID simulated annealing results for the solar 
equator. In Panel (a) The crosses are the observed splitting com- 
binations and the circles are those obtained by the fit. Panel (b) 
shows the normalised residuals. 

scatter, which is reflected in comparatively large errors even 
though the errors in splitting coefficients is lowest for these 
latitudes. This could be due to the fact that the form of trend 
assumed in the tachocline model defined by equation (17) 
is not sufficient to model the variation in rotation rate. As 
can be seen from the contour diagram in Fig. 3, there is 
some variation at low latitudes in the convection zone, which 
perhaps cannot be modelled properly by a linear trend. 

The annealing results also indicate that the thickness of 
the tachocline is very small at low latitudes. The thickness 
appears to increase at higher latitudes though the estimated 
errors also increase and it is not clear if the increase is sig- 
nificant. Similarly, it is not clear if shift in the tachocline 
position with latitude is also significant, though in general 
once again the tachocline appears to shift to slightly larger 
radial distance at higher latitudes. 

We feel that there should be yet another independent 
test of the significance of variation in tachocline position and 
thickness with latitude. We have therefore tried a 2D anneal- 
ing fit to simultaneously determine all the parameters at all 
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Figure 11. The fit to the first three splittings coefficients and 
the normalized residuals of the fits obtained by the 2D simulated 
annealing method assuming that there is no latitudinal variation 
in the position and thickness of the tachocline. 

latitudes. Since this fit involves 12 parameters as defined by 
equations (17), (18), it is not clear if all of them are required. 
For verifying this we start with the simplest situation where 
the position and width of tachocline are independent of lati- 
tude (i.e., rd3 = 0, W3 = 0). This fit yields the mean position 
of tachocline as O.7O38P0 and a half-width of 0.0187P© and 
the fit is shown in Fig. 11. These values are consistent with 
the results obtained by Basu (1997), though the half- width 
is slightly higher than the value found by Basu (1997). Some 
of the difference may arise because while Basu (1997) used 
only the splitting coefficient C3 to determine the tachocline 
parameters, in this work we are using the first three splitting 
coefficients. If this difference is significant it may imply a lat- 
itudinal variation in tachocline position or width. Further, 
this fit has a x 2 — 1-108 per degree of freedom and the fit 
appears to be reasonably good. When all 12 parameters are 
included in the fit the minimum \ 2 P er degree of freedom 
reduces to 1.088. Thus the reduction in \ 2 is only marginal 
and it is tempting to conclude that the GONG months 4-14 
data are consistent with no latitudinal variation in position 
or width of the tachocline. Of course, we can not rule out 
a small variation in the tachocline properties with latitude. 
Monte-Carlo simulation with 12 parameters yields the fol- 
lowing results for the tachocline: 

r d = [(0.6991 ± 0.0099) + (0.0030 ± 0.0061)P 3 (#)]#©, 
w = [(0.0084 ± 0.0072) + (0.0047 ± O.OO42)P 3 (6O]P , 
6Q = (-1.83 ± 2.18) - (22.71 ± l.Ol)P 3 (0) 
- (3.88±0.45)P 5 (6>) nHz, 
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It is clear that the variation in position is at a level of (l/2)cr 
only, while the thickness at all latitudes is comparable to the 
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Figure 12. Same as Fig. 11, but with the possible latitudinal 
variation in the position and thickness of the tachocline taken 
into account 

error estimates and as such it is not clear if the variation in 
thickness is significant. Similarly, the first component of 5Q 
is also not significant. In fact, we find that if this parameter 
is set to zero, the fit is more stable in the sense that the 
X 2 reduces to acceptable level in most of the attempts with 
simulated annealing and as such we prefer to use that fit. 
The minimum value of x P er degree of freedom in this case 
is 1.105 and the resulting parameters of the tachocline are 

r d = (0.6947 + 0.0035P3(6»))P Q , 

w = (0.0067 + 0.0014P 3 (6»))P Q , (22) 
= -21.11P 3 (0) - 2.96P 5 (0) nHz, 

while the fit is shown in Fig. 12. We adopt these values for 
the tachocline model for use in the next section and for com- 
parison with other estimates. It is clear that addition of two 
parameters determining the variation in properties with lat- 
itude does not improve the fit significantly and hence there 
is no compelling reason to believe that there is any vari- 
ation in position or thickness of tachocline with latitude. 
Nevertheless, all the results which attempt to determine the 
variation find an increase in thickness with latitude and a 
shift outwards in the mean position of tachocline with in- 
creasing latitude. Although this variation does not appear to 
be significant in terms of the expected errors, a small varia- 
tion in properties of the tachocline with latitude cannot be 
ruled out. 

Fig. 13 summarizes the results from all the techniques. 
This figure shows the variation in jump, position and thick- 
ness of the tachocline as obtained by different techniques 
from the GONG months 4-14 data. We have chosen the re- 
sults from calibration method using the models with trend 
for this purpose. It is clear that there is a general agreement 
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Figure 13. A summary of the tachocline results obtained using 
the GONG months 4—14 data. Panels (a),(b) and (c) show the 
results of the jump, position and thickness respectively. In each 
panel the continuous line is the results of the 2D annealing (equa- 
tion 18) with the la error bounds shown as the dotted lines. The 
circles show the results of the calibration method and the trian- 
gles are the results obtained by ID annealing. The symbols are 
displaced by ±0.5° about the true latitude for the sake of clar- 
ity. The dashed line in panels (b) and (c) mark the mean values 
found by Basu (1997), while the lcr error bounds are shown as 
dot-dashed lines. 

between the three independent results at most latitudes and 
all these results point to an increase in thickness with lat- 
itude as well as an outward radial shift in the position of 
tachocline with latitude. This figure also shows the mean po- 
sition and width of tachocline as determined by Basu (1997) 
and it appears that all the results are also consistent with 
no latitudinal variation in position or width of tachocline. 
Clearly, better data is required to find any possible variation 
in tachocline properties with latitude. 



4 INVERSION AFTER REMOVING THE 
TACHOCLINE SIGNAL 

The smoothing used in our inversion procedure tends to 
smooth out the steep variation in rotation rate in the 
tachocline. Apart from this it may also introduce some rip- 
ples away from the position of tachocline, a feature that is 
reminiscent of the Gibbs phenomenon in Fourier transform. 
In order to overcome this limitation and to use the tachocline 
parameters determined in the earlier section for improving 
the results of inversion, we attempt an inversion after remov- 
ing the signal from tachocline. For this purpose we adopt 
tachocline parameters given by equation (22), leaving aside 
the smooth part and perform a forward calculation to ob- 
tain the splitting coefficients for this model. These splittings 




Figure 14. A contour diagram of the solar rotation rate as ob- 
tained by 1.5D inversion of the GONG months 4-14 data after 
removal of the tachocline. The format is the same as that for 
Fig. 3. 

are then subtracted from the observed splittings before in- 
verting the data. The rotation rate in the tachocline model 
is then added to the inverted profile to obtain the actual 
rotation rate in the Sun. 

The results in Fig. 14 display the contour diagram for 
the rotation rate obtained using 1.5D inversion, while Fig. 15 
shows the rotation rate as a function of radial distance at se- 
lected latitudes for the inversions with and without removal 
of tachocline. From this figure it appears that at the 30° 
latitude the results obtained using the 2D inversion after re- 
moving the tachocline are different from other results. The 
reason for this difference are not altogether clear, though 
it appears to manifest at latitudes around the region where 
the radial gradient in the tachocline changes sign. The re- 
sults obtained after removing the tachocline appear to be 
smoother in general, which is of course, not surprising. Also, 
the fact that the simulated annealing was able to obtain a 
good fit over a large fraction of radial distance using only 
12 parameters suggests that there is probably not much 
structure in the rotation rate as a function of latitude or 
radial distance (apart from tachocline itself) in the region 
0.6 < r / Rq < 0.9. However, from the 2D annealing fits 
there are reasons to believe that these 12 parameters are 
not completely adequate to represent the trends at all lati- 
tudes adequately. Nevertheless, it is likely that some features 
seen in the lower part of CZ in certain inversion results (e.g., 
Figs. 3,4) are due to the influence of tachocline magnified by 
errors in data. 



5 CONCLUSIONS 

In this work we have attempted to infer the rotation rate 
inside the Sun using observed p-mode frequency splittings. 
We have used a 1.5D inversion technique to invert the data 
in the form of splitting coefficients based on a regularized 
least squares method with iterative refinement. We have in- 
vestigated the influence of using different prescription for 
smoothing to find that the differences are comparable to ex- 
pected errors in inversion results. Similarly, an adoption of 
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Figure 15. A comparison of the inversion results with and with- 
out removal of tachoclinc. Continuous line shows the results ob- 
tained using 1.5D inversion after removing the tachocline contri- 
bution, the dotted line is the 1.5D result of normal inversion. Simi- 
larly, the dashed line show the results of 2D inversion of individual 
splittings after removing the tachoclinc while the dot-dashed line 
shows the results from normal 2D inversion. 

different sets of observed splitting coefficients shows that the 
difference in various results is again roughly consistent with 
error estimates based on quoted errors in observed split- 
tings. Apart from the 1.5D inversion technique, we have also 
attempted 2D inversion for both the individual splittings 
D ny t,m and the splitting coefficients ci n ' e \ All these results 
agree reasonably well inside the convection zone, although 
in the deep interior there are some differences between dif- 
ferent inversion results. We believe these reflect our inability 
to obtain reliable estimate of rotation rate in the core due to 
possible systematic errors in splittings for low degree modes. 

It is clear from our results that the surface differential 
rotation persists through the solar convection zone, while 
below the base of convection zone the rotation rate ap- 
pears to be relatively independent of latitude. The tran- 
sition around the base of the convection zone may not be 
resolved by the inversion results. The core appears to be ro- 
tating slower than the surface equatorial rotation rate as also 
found by Tomczyk, Schou & Thompson (1995) and Elsworth 
et al. (1995). The rotation rate in solar core has been re- 
cently discussed by Gizon et al. (1997) and Rabello-Soares 
et al. (1997) and it appears that there is significant variation 
in the inferred rotation rate in the solar core from different 
helioseismic data. The constraints on the rotation rate in 
the core in the form of the boundary conditions (8) , in the 
1.5D inversion make the core rotate essentially as a rigid 
body and reduces the discrepancy in rotation rate inferred 
using different data sets. This may be artificial, but unfortu- 
nately we do not have enough data to resolve any variation 



in rotation rate in the core and this is probably the simplest 
assumption that we can make (Chaplin et al. 1996). Even 
a small error in splittings for low degree modes can change 
the inverted rotation rate in the core significantly if no con- 
straints are applied. Further, as can be seen from Fig. 1, the 
results outside the core are not too sensitive to the point 
where these boundary conditions are applied or the form of 
smoothing used. The main problem with reliable inversion 
of rotation rate in the core is that the low degree modes that 
penetrate into the core have large errors in the splitting co- 
efficients (possibly including some systematic errors), and 
this allows for a wide range of rotation profiles in the core. 
It is clear that better quality data are required to make any 
reliable estimate of rotation rate in the solar core. 

Further, there is a distinct shear layer just underneath 
the solar surface where the rotation rate increases with 
depth. There is also a hint of this shear layer becoming 
less pronounced with latitude, in the sense that the radial 
variation of the rotation rate tends to diminish with in- 
creasing latitude. The rotation rate has a maximum around 
r = 0.95.R© at most latitudes, except possibly close to the 
poles. This feature can also be seen in the raw splitting 
data (Fig. 7). The inferred rotation rate at the solar sur- 
face agrees reasonably well with that obtained from Doppler 
measurements (Snodgrass 1992). Some differences between 
the Doppler measurement and inverted profile (Fig. 6) at 
higher latitudes could be because the Doppler measurements 
determine only terms up to cos 4 6 in expansion of rotation 
rate. The higher order terms included in inversions will con- 
tribute significantly at high latitudes. 

The tachocline or the shear layer near the base of con- 
vection zone has been studied in detail using forward tech- 
niques to ascertain possible latitudinal variation in its prop- 
erties. Since the tachocline cannot be resolved by inversions, 
we have used a number of forward modelling techniques to 
study the tachocline. With the present data, we do not find 
any compelling evidence for any variation in the position or 
thickness of the tachocline with latitude. The mean position 
and thickness of tachocline is found to be consistent with the 
values found by Basu (1997). Our results suggest that the 
thickness increases marginally with latitude and the location 
of tachocline also appears to shift outwards with increasing 
latitude, but the difference between the position and thick- 
ness of the tachocline at the solar equator and at a latitude 
of 60° is comparable to the estimated errors. It is therefore, 
not clear if the variation is really significant. Similarly, the 
2D annealing results also show that the variations in the po- 
sition and thickness are less than the respective error bars. 
Further, there is no significant reduction in the \ 2 f° r the 
fit when the variation of position and width with latitude is 
included in the 2D fits. Taking the errors into account we 
believe that we get an upper limit of 0.03J?q for the varia- 
tion in the position of the tachocline and about O.O27? for 
the variation of the thickness. Clearly, more precise data are 
needed to make a better estimate of the parameters reliably. 

From our results it appears that the tachocline is cen- 
tered at a depth which is below the base of the solar convec- 
tion zone at all latitudes. If the latitudinal variation shown 
by our results is real then at low latitudes most of the varia- 
tion in rotation rate within the tachocline occurs below the 
CZ base, while at high latitudes part of the tachocline may 
extend into the convection zone. 
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